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Abstract. We analyze the role of the interplay between on-site interaction and 
inhomogeneous diffusion on the phenomenon of condensation in the zero-range 
process. We predict a universal phase diagram in the plane of two exponents, 
respectively characterizing the interactions and the diffusion disorder. The most 
prominent outcome is the existence of an extended condensed phase. In the 
latter phase, which originates as a result of the combined effects of strong enough 
interaction and weak enough disorder, a typical high-density configuration has a 
unique condensate on top of a critical background, but the condensate may be 
located at any site of a large hosting set of favored sites, whose size grows sub- 
extensively. The novel extended condensed phase thus interpolates continuously 
between the two scenarios associated so far with the condensation transition, 
namely spontaneous symmetry breaking and explicit symmetry breaking. 



E-mail: Claude . godrecheOcea. fr, jean-marc . Iuck0cea.fr 



Condensation in the ZRP: an interplay between interaction and disorder 



2 



1. Introduction 

The essence of the zero-range process (ZRP) can be best understood by returning 
to the original work of Spitzer pQ. The starting point is to consider an assembly of 
independent particles diffusing on a graph, i.e., a finite connected set of M sites. The 
particles, to which we can also refer to as random walkers, are also in finite number, N. 
Choosing a continuous-time description of the process, any individual particle hops 
from site m to site n with a rate w mn . In the course of time, the probabilities for 
the particle to be at the various sites m relax to the (unique) stationary occupation 
probabilities q m , obeying the master equation 

M M 

} J w nm q„ =q m /, w mn . (1.1) 

n— 1 n— 1 

For the assembly of independent random walkers, if at time t the occupation of site m 
is equal to k, then, of course, the chance of seeing a particle hop from this site to site n 
is enhanced by a factor fc, which amounts to multiplying the hopping rate w rnn by k. 

The intention stated by Spitzer was to depart from the field of independent 
random walkers by introducing a form of interaction between them. In Spitzer's 
original definition of the zero-range process, the rate k w mn for the hopping of a 
particle from site m, with occupation equal to k, to site n, which holds if the particles 
are independent, is changed to the product 

W(m,n,k) = w mn u k , (1.2) 

where the factor Uk only depends on the occupation k of the departure site m at the 
current time. This dependence gives an explicit content to the denomination of the 
process as 'zero-range'. This very definition of the ZRP demonstrates the importance 
and respective roles of the two facets of this process: diffusion on the one hand, 
interaction on the other hand. The first one is encoded in the single-particle hopping 
rate w mn , while the second one is encoded in the on-site interaction factor u/~. Either 
factor of the product (jl.2l) can be separately normalized in an arbitrary way (see 
below (|2.5p ). This justifies that either factor will sometimes be referred to as a 'rate'. 

The interest of the statistical physics community for the ZRP came long after the 
publication of [1] and the subsequent works of mathematicians on the subject [21 E]. 
The popularity of the ZRP was triggered by the ubiquity of the phenomenon of 
condensation in a number of processes related in one way or another to the ZRP (for 
reviews, see [4] El E]). The condensation transition manifests itself by the existence of 
a critical density p c , with the following meaning. Consider the stationary state of the 
model in the thermodynamic limit where M and N ~ > oo, at fixed density p — N/M. 
If the density is high enough {p > p c ) 1 a macroscopic (extensive) number of excess 
particles, of order 

A = JV — Mp c = M(p — p c ), (1.3) 

concentrates on a single site with very high probability. 

In the present work we analyze the role of the interplay between the two facets 
of the ZRP, diffusion and interaction, on the phenomenon of condensation. More 
precisely, we consider the following model, defined on a finite system consisting of N 
particles living on M sites. 

• We characterize the diffusion disorder, i.e., the inhomogeneity of the diffu- 
sion process, by modeling the single-particle weights (stationary occupation 
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probabilities) q rn as i.i.d. random variables drawn from a continuous distri- 
bution f(q), chosen as 

/(<?) = c(l - qf- 1 (0 < g < 1), (1.4) 

with a power-law singularity as q — > 1 with an arbitrary exponent c > 0, to be 
referred to as the disorder exponent. The shape of the distribution /(g) in the 
vicinity of the maximum (supposed to be finite and set here to q m &x = 1) will 
play a central role in the analysis of the condensation transition. 

We choose for the factor Uk, which encodes the interaction between particles 
sitting at the same site, the form 

b , \ 

u k = l + r, (1.5) 

where & is a parameter describing the strength of on-site interactions, to be 
referred to as the interaction exponent. 



The above model will be further motivated in section 12.41 It contains as limiting 
cases two known models exhibiting a condensation phenomenon (see references below) . 

(i) Homogeneous ZRP. If diffusion is homogeneous (q m — 1, i.e., f(q) = 5{q — 1), 
or formally c = 0), the process with rate (II. 5j) is known to have a condensation 
transition at the critical density 

Pc - ^- 2 , (i.6) 

whenever the interactions arc sufficiently attractive (b > 2). In this situation, the 
condensate can sit on any site with equal probabilities in the stationary state of 
the system. This condensation scenario therefore corresponds to a spontaneous 
symmetry breaking (SSB) of the symmetry between the sites. 

(ii) Occupation-independent inhomogeneous ZRP. If the rate uu is occupation- 
independent (uk = 1, i.e., 6 = 0), the process with inhomogeneous diffusion 
described by the distribution (|1.4j) is also known to have a condensation transition, 
at the critical density 

Pc = ^ (1.7) 

whenever the disorder is sufficiently strong (c > 1). In this situation the 
condensate is typically located on the site with the largest single-particle weight. 
This corresponds to a scenario of explicit symmetry breaking (ESB). 

The full model defined by both diffusion disorder (|1.4[) and interaction (|1.5|) 
interpolates between the two limiting cases recalled above. The aim of this work 
is to investigate to what extent the condensation transition of the homogeneous ZRP 
defined in (i) is altered by the presence of inhomogeneous diffusion, i.e., of the non- 
trivial distribution (|1.4|) of the single-particle weights q m , or, conversely, how the 
condensation transition of the inhomogeneous ZRP defined in (ii) is altered by the 
presence of an interaction corresponding to the rate (jl.5l) . 

An informal account of our results is as follows. We show in figure Q] the phase 
diagram of our model in the plane of the exponents b (characterizing interaction) and c 
(characterizing disorder). In the fluid phase, there is no condensate at any density. 
In the localized condensed phase, for any p > p c , a condensate of size A (see (|1.3|0 
lives on the site with the largest single-particle weight with very high probability. In 
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the extended condensed phase, for any p > p c , the condensate is still unique with very 
high probability, but it may be located at any site of a large hosting set of favored 
sites, whose size R grows sub-extensively as 

R~£~j^. (1.8) 

A c {p-Pc) c 

The existence of an extended condensed phase is the most salient novel feature put 
forward in the present work. The latter phase implements a continuous interpolation 
between the two symmetry breaking scenarios recalled above, namely SSB in the 
homogeneous ZRP (c = 0, b > 2), and ESB in the occupation- independent 
inhomogeneous ZRP (6 = 0, c > 1). These two limiting scenarios are consistently 
recovered as R = M and R = 1, respectively. The ESB scheme more generally applies 
to the ZRP with c > 1 and any b. Figure [2] shows a more refined phase diagram, where 
the condensed phases are further divided into a regime of normal fluctuations, where 
the mean square occupation /i c at criticality is finite, so that the size fluctuations of 
the condensate around its mean size A are Gaussian and grow as M 1 ' 2 , and a regime 
of anomalous fluctuations, where p, c is divergent, so that size fluctuations are not 
Gaussian and grow faster than M 1 / 2 . The phase diagrams depicted in figures [1] and [2] 
are universal, in the sense that they hold for any rate whose asymptotic behavior at 
large k is of the form w 1 + b/k, and for any distribution of single-particle weights 
exhibiting a finite maximum g max and a power- law singularity f(q) ~ (cjWx — <?) c l 
near this maximum. 



Localized condensate 



Fluid 




Figure 1. Universal phase diagram of the inhomogeneous ZRP in the b-c 
plane, with its fluid phase and its two condensed phases: a localized one 
(c > 1) and a novel extended one (c < 1). 



Before proceeding, let us review the relevant literature. 

The homogeneous ZRP with rate (|1.5j) has been studied in a long series of works, 
and is by now well understood [H [5J [6] . The main features of the statics of the model 
without reference to a dynamical process were first analyzed in [?Jj] . The underlying 
dynamical process was introduced in [8], where the coarsening dynamics leading to 

| The equilibrium state of the static model considered in [7] is actually slightly different from the 
stationary state of the homogeneous ZRP with rate l|1.5| l ■ However the universal static properties of 
the two models are the same. The same holds for the dynamics, i.e., the universal properties of the 
dynamics of the homogeneous ZRP with rate 11.51 1 or of the process considered in [8] [9] are the same. 
This was underlined in |10| . 
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Figure 2. More detailed phase diagram, showing the internal structure of 
the condensed phases. The size fluctuations of the condensate around its 
mean A are either normal (finite mean square occupation), or anomalous 
(divergent mean square occupation). 



condensation was studied in the mean-field geometry of the complete graph. This 
study was subsequently pursued in [5j. The relationship of the equilibrium state of 
the model described in [7] to the stationary state of a ZRP was recognized in [IT] 
and a first analysis of the criterion on the form of the rate Uk for condensation to 
occur was given. This analysis was subsequently completed in J3J. Independently, a 
similar analysis was performed in [12], devoted to a study of the size of the condensate. 
Further studies involve the analysis of the coarsening process leading to condensation 
in the one-dimensional case [TO] D~3] , or on the complete graph [10] , the proof of the 
uniqueness of the condensate and the equivalence of ensembles [13], the analysis of 
the motion of the condensate at stationarity [T4] , the analysis of finite-size effects in 
the canonical ensemble [15] , and the relationship of the condensation phenomenon to 
extreme- value statistics [IHj. More mathematical works on similar subjects can be 

found in [13 HEIDI. 

The occupation-independent ZRP with rate Uk = 1 and inhomogeneous diffusion 
described by (|1.4j) can be traced back to studies of one-dimensional disordered totally 
asymmetric exclusion models [20] [21] [22] , where the connections with the ZRP [20] [21] 
and the relationship to Bosc-Einstcin condensation |22j were recognized. Subsequently 
the dynamics of the same model was analyzed in [23] , both at stationary and during 
the process of the formation of the condensate. 

The references cited so far led to the current understanding of the two limiting 
prototypical cases of the ZRP referred to above. We now turn to other works 
which deal in some way or another with the interplay of interaction and disorder 
in the ZRP. The simplest situation of an inhomogeneous ZRP is that of an attractive 
impurity in a homogeneous background: one site is singled out by its weight q > 1, 
while all the other sites keep q = 1. In this situation condensation occurs even 
for the occupation-independent ZRP (uk = 1) [H EI]- The situation where the 
diffusion disorder originates in the structure of the underlying graph have also been 
analyzed [521 ESI EI], the main emphasis being on complex (scale-free) networks. 
In [55] complete condensation, i.e., p c = 0, has been shown to generally occur 
on such networks. References [26] also deal with the statics and the dynamics 
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of the occupation-independent ZRP on complex networks. The phenomenon of 
condensation in disordered urn models has been explored in [55]. The more general 
case of an inhomogeneous ZRP where disorder enters the interaction part of the rate 
itself has been recently considered in [29] . Finally, the occupation-independent one- 
dimensional ZRP with partially asymmetric disordered diffusion has been analyzed 
in [30) . Reference [31] presents a brief review of the more general area of driven 
diffusive systems with disorder. 

The setup of the present paper is as follows. In sectional we present some general 
formalism and emphasize two special situations, namely Uk = k (independent random 
walkers) and uu = 1 (the occupation-independent case, which can be mapped onto an 
ideal Bose gas). We then review, in section[3l the two limiting cases of interest recalled 
above, namely the homogeneous ZRP with rate (ll.5[) and the occupation-independent 
inhomogeneous ZRP. Section 2] contains our theoretical investigations, whose main 
predictions are summarized in the universal phase diagrams depicted in figures [T] 
and[2J Our quantitative predictions are then illustrated by means of numerical results 
in section [5j Section [6] presents an application of our findings to the case of complex 
networks. Section [7] contains a brief discussion, whereas more technical aspects are 
exposed in three appendices: |Appcndix A| contains an investigation of the ZRP on 
two sites, while |Appcndix B| is devoted to the universal fluctuations of the hosting 
probabilities in the borderline situation (c = I), and |Appcndix Cj to an effective model 
describing the relevant rare events in the localized regime (c > I). 



2. General formalism 

2.1. Stationary state 

We start by a reminder of the definition of the ZRP and of the properties of its 
stationary state. Consider a finite connected graph of M sites (labelled m = 1, . . . , M), 
on which particles hop from site to site in continuous time. We denote by N m the 
instantaneous occupation of site to, i.e., the number of particles living at site m at 
the current time. The total number of particles in the system, 

M 

N=J2 N m, (2.1) 

m— 1 

is conserved by the dynamics. 

According to Spitzer's original work pQ, a ZRP is defined by the following 
expression (see (|1.2[l ) 

W(m, n, k) = w mn u k (2.2) 

for the rate for a particle to hop from site m (with occupation N m — k) to site n, 
where w mn describes the diffusion of a single particle on the graph, and u k encodes 
the interaction between particles sitting at the same site. In the present work we focus 
our attention on an inhomogeneous ZRP defined by (|2.2[) . i.e., with diffusion disorder. 
For completeness, we mention the following more general definition of a ZRP [5], 

W(m, U, N m ) = WmnUm.k, (2.3) 

where the factor u m ^ encoding the interaction between the k particles sitting at site m 
depends explicitly on both m and k. References [29] contain an investigation of such 
a situation, where quenched disorder affects the interactions between particles, and 
not only the diffusion. 
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The fundamental property of the ZRP is the simple structure of its stationary 
state: the distribution of the occupations is factorized. Such a distribution is 
referred to as a product measure. The particles being considered as indistinguishable, 
a configuration of the system is entirely defined by the occupations {N m }. The 
stationary probability P({N m }) of any configuration is given by a product of 
elementary weights, associated with the occupations of each site. For the hopping 
rate (jT2]l (or ([2~2]) ). we have 

m / m \ 

P({N m }) = — J] q%» PNm S( N mi N . (2.4) 

' m— 1 \m— 1 / 

The Kronecker delta function ensures the condition (|2.ip . whereas the normalisation 
factor Zm,n (see (|2.6p ) plays the role of a partition function, and the weight associated 
with each site m is itself a product of the two factors, q^ m and PN m ■ The latter factor 
is given by 

Po = l, Pk = • (2.5) 

Ul . . .life 



As long as one is only interested in static properties, either factor of the rate (( 1 . 2|) 
(or (|2.2I) ) can be normalized in an arbitrary way. Indeed, if all the hopping rates w mn 
arc multiplied by the same constant A, the time scale of the dynamical process 
is changed as t — > t/A, but the single-particle weights q m obeying the stationary 
equation are left unchanged. Similarly, if all the factors u k are multiplied by 

the same constant B, p k is changed to pk/B k , the partition function Zm.n gets 
multiplied by 1/B N , and so the result (|2.4[) is unchanged. Furthermore, the single- 
particle weights themselves can also be normalized in an arbitrary way. Indeed, if all 
the q m are multiplied by the same constant C, Zm,n gets multiplied by C , and so 
the result (I2.4[) is again unchanged. 

Canonical ensemble. In the canonical ensemble, the total number N of particles 
is fixed. The stationary probability P({N m }) of any configuration is given by the 
expression (|2.4p above, where the partition function 



M I M \ 

\N m }m=l \m=l ) 



Zm.n = E 11 C m PN m 5[ ^ N m ,N ) (2.6) 

{N m } m=l 

depends on all the q m and all the pk- 

Singling out the sum over N m = k in (12.61) leads to the recursion relation 

N 

Zm.n = ^q } mPk Z M -i,N-k(4 m ) (2-7) 

between partition functions of systems with successive sizes. The argument $ m means 
that ZM-i,N-k{i m ) depend on the (M — 1) single-particle weights q n for n ^ m. 
Likewise the distribution of the occupation of site m is obtained from (|2.4p as 

t o kfAr M k Z M -l,N-k{i m ) , n Q , 

f m ,k = Prob{iV m = k\ = q m p k —. (2.8) 

AM, AT 

Hence the mean occupation (local density) at site m reads 

Pm = (N m ) = kf rn k = gm - M ' N = q mT ^—lnZ M .N. (2.9) 
^ Z M ,N oq m dq m 
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The partition function Zm,n being a homogeneous function of degree N of all the q m , 
the sum rule 

M 

Y,Pm=N (2.10) 
m— 1 

is ensured by Euler's identity for homogeneous functions: 

E ftn = NZm,n- (2.11) 

m—1 

By using a contour-integral representation of the Kronecker delta function, we 
can recast the partition function (12.61) into the form 



■>M,N 



2mz N + 

where the grand partition function reads 



T Zm(z), (2.12) 



M 



z M {z) = x: z M , N z N = n p(^ m ), (2.i3) 

N>0 ro=l 

and 

fc>0 

Grand- canonical ensemble. The representation (|2 . 12|) of the partition function sug- 
gests to alternatively consider the grand-canonical ensemble, where the fugacity z 
is fixed, whereas the total number N of particles in the system fluctuates. In that 
ensemble the stationary distribution of the occupations is strictly a product measure: 

M M 

p G °({N m }) = -—r n ( Z q m ) N - PNm = n (2.i5) 

m—1 m—1 

where the normalized factor f^ k gives the distribution of the occupation of site m: 
f^ k = Prob{^ c = k} = (2.16) 



Thus the mean occupation at site m reads 

P m C = (N* C ) = £ = Zq "£' {Zq : i) = ^ ^P(.zq m ), (2.17) 

where the accent denotes differentiation, and finally the mean density of the system is 



r( . (N GC ) 1 rr zZ'(z) z d , . . 

gc = \ / = _ gc = m}_) = _ z 2 18 

M M^ Hm MZm(z) Mdz MK ' V ' 



2.2. The case of independent random walkers (uk = k). 

We now turn to the simplest example of a ZRP, corresponding to Uk = k. As recalled 
in the introduction, this situation describes independent random walkers. We have 

Pk = i P(z) = e z . (2.19) 
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Grand- canonical ensemble. The grand partition function reads 

Z M {z)=e Mz ~\ (2.20) 
where we have introduced the notation 
1 M 

«=mE^ ( 2 - 21 ) 

m— 1 

for the mean one- particle weight. The occupations therefore follow the Poissonian law 

/£S = e-*»^|f£, (2.22) 
and the local densities read 

pl G = zq m . (2.23) 
Thus the mean density and the fugacity are related through 

p GC = zq. (2.24) 

As a consequence, for any density p, there exists a value of the fugacity, namely 
z = p/q, such that p GC (z) = p. The system can sustain any density of particles and 
is therefore always in a fluid phase. 



Canonical ensemble. The canonical partition function is 

{Mq) N 



(2.25) 
(2.26) 



* M > N m 

and the local densities read 

_ Nqrn 
Pm ~ Mq ■ 

The canonical joint distribution of the occupations is multinomial: 

1 /V! M 

p ««)=(M^v^r— IT^ m ' ( 2 - 27 ) 

n m_i 

m— 1 

as it should be, since independent random walkers behave as independent classical 
particles. 

2.3. The occupation-independent case (uk = lj. 

The occupation-independent ZRP with Uk — 1 is another simple example, yet richer 
because the random walkers are no longer independent. We have then 

Pk = 1, P(z) = (2.28) 
l — z 
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Grand- canonical ensemble. The grand partition function reads 

M 



= II T^TT- ( 2 - 29 ) 

m—1 

The occupations therefore follow the geometric law 

= (1 - zq m )(zq m ) k , (2.30) 
and the local densities read 

GC _ ZC bn 



r m 1 - zq m 

Thus the mean density and the fugacity are related through 

M 



(2.31) 



m—1 

For a homogeneous system (g m = 1), th above relation simplifies to 

p GC = t4^' ( 2 - 33 ) 

which has a solution z(p ) for any value of the density, and hence the system is again 
in a fluid phase. However, if the system is inhomogeneous, local densities are given 
by (I2.31|) . and a condensation transition may occur. The simplest case is that of an 
attractive impurity [4j [24] , where one of the q m is larger than the common value of 
all the other ones. A condensation transition also occurs in the thermodynamic limit 
when the q m are random and distributed according to (|1.4[) with c > 1. The latter 
phenomenon is addressed in section 13.21 

It is worth emphasizing the analogy between the statics of the occupation- 
independent ZRP and that of an ideal Bose gas [22]. Setting z — e^, where j3 
is the inverse temperature and /i the chemical potential, and q m — e~ l3Em , the 
expression (|2.31j) for the local particle density at site m identifies with the Bose 
occupation factor 

/Bose(-Em) = ^E^JT) _ j ( 2 - 34 ) 

of a fictitious quantum-mechanical level with energy E m . Identifying the connection 
between classical stochastic processes and the ideal Bose gas dates back to [32] . 

The distribution (|1.4p of the one-particle weights translates to a distribution of 
energy levels (density of states) D(E) vanishing as D(E) ~ E c ^ 1 as the energy E > 
approaches zero. The latter power-law singularity is to be put in perspective with the 
density of states of a free particle in a box in d-dimensional space, i.e., D(E) - E d l 2 - 1 . 
This parallel leads to the identification c = d/2 between the disorder exponent c and 
half the dimensionality of the space where the Bose gas lives. The condition c > 1 
for having a condensation transition [2U [22], recalled above and to be worked out 
in section 13. 2[ is therefore equivalent to the well-known result that an ideal Bose gas 
exhibits a Bosc-Einstein condensation in dimension d > 2 only. 
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Canonical ensemble. The canonical partition function can be evaluated from the 
representation (12.12)) as 

M 



z 



M,N = 



2mz N + 1 11 1 - zq. 

m—l 



Evaluating the contour integral by the method of residues, we get 

M 

Zm,N — Q.mQmi (2.36) 



with 



n 



Qm Qn 

n^rn 



(2.37) 



The local densities can then be obtained by differentiating the above formulas 
according to (|2.9p . We thus obtain 

Pm = ^— \ NqZQrn + V ^^ + ^Qn \ ^ 

^m,n \ '-r q n — q-m I 

\ n^m J 

The expressions (|2.36l) and (|2 .38[) are much more intricate than their counter- 
parts (|2.25j) and (|2.26|) in the case of independent random walkers. This difference 
can be attributed to the fact that the canonical ensemble is the natural framework 
to describe identical classical particles (hence the simplicity of the multinomial joint 
distribution (|2.27|) ). whereas the grand-canonical ensemble is the only natural one 
when dealing with quantum-mechanical energy levels (hence the analogy between 
the expression (|2.31[) for the grand-canonical local density and the Bose occupation 
factor (|2~34| ). 



2.4- On the choice of single-particle weights made in this work 

Let us come back to the specific model considered in the present work. For a given 
ZRP, defined by the rates Uk and w mn , the stationary properties only rely on the 
knowledge of the weights pk and q m . The former ones are easily deduced from the 
rates Uk (see (|2.5j) ) . In contrast, the stationary master equation (|1.1[) cannot be solved 
in closed form in general. The problem indeed belongs to the rich and complex area 
of diffusion in random media [33 . In this work we have adopted the alternative 
viewpoint of modeling diffusion disorder by considering the single-particle weights q m 
as i.i.d. random variables drawn from the distribution (|1.4|) . 

Let us take the example of the mean-field geometry of a completely connected 
graph over M sites to discuss the issue. Consider the inhomogeneous diffusion where 
the particle leaves site m with some site-dependent rate A m , and jumps to a random 
arrival site n, chosen uniformly over the system. We have therefore w mn — X m /M. It 
can be checked that the stationary probabilities read q m = 1/A TO , up to normalization. 
In the more general situation of separable rates of the form w rnn — A m /i n , we have 
qm = A* m /A TO , again up to normalization. The stationary state thus obtained is 
an equilibrium state, obeying the condition of detailed balance. The mean current 
from site m to site n, J mn = w mn q m = [i m Hn, up to normalization, is indeed 
symmetric. Another example is the case of the totally asymmetric diffusion on a one- 
dimensional lattice, either infinite of finite with periodic boundary conditions, such 
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that w mn = \ m $n,m+i- Here again, we have q m = 1 / A m , up to normalization. The 
above examples of inhomogeneous diffusion in the mean-field geometry and in the one- 
dimensional totally asymmetric case demonstrate that the single-particle weights q m 
can take any prescribed values. This somehow justifies our choice of modeling them by 
i.i.d. random variables. The rationale for the choice (jl.4l) of the distribution f(q) of 
single-particle weights is as follows. The properties of the condensate are entirely 
dictated by the behavior of this distribution f(q) near its maximum g max . The 
situation where g max is infinite is somehow pathological (section [5] contains a study 
of a physically motivated case of this type, namely complex networks). In the case 
where q ma x is finite (and can then be set equal to unity) , the situation of a power-law 
singularity turns out to be the appropriate framework to disclose the phase diagram. 
Finally, the pure power law (jl.4l) has been chosen merely for simplicity. 



3. Two limiting case studies 

As mentioned in the introduction, the full model defined by (|1.4[1 and (| 1 . 5|) contains 
two limiting prototypical cases which are known to exhibit a condensation transition. 
We recall here the properties of these models which will be relevant for the study of 
the full model, to be presented in section 2) 



3.1. Homogeneous ZRP: spontaneous symmetry breaking 

The homogeneous ZRP (see H |] for reviews) corresponds to the case where 
the single-particle weights q m are uniform, i.e., do not depend on the site. This 
homogeneity property holds for the complete graph with uniform hopping rates, and 
more generally whenever all sites are symmetry-related and thus equivalent. This is 
automatically satisfied if the sites are related to each other by translation invariance, 
like e.g. for usual symmetric or biased random walk on finite lattices with periodic 
boundary conditions. 

Consider the homogeneous ZRP with rate (|1.5I) . Setting q m = 1, we have 

Z M {z)=P{z) M . (3.1) 

As a consequence, in the grand-canonical ensemble, (|2.16[) and (|2.17l) respectively 
become 

k 

h -Pjz) {S - 2) 

and 

" G ° = W < 3 - 3 > 

In the canonical ensemble, in the thermodynamic limit (M — > oo, N — > oo, p = 
N/M fixed) , the condition (|3.3I) is recovered by estimating the contour integral (I2.12[) 
by the saddle-point method. This shows that both statistical ensembles are equivalent 
in the thermodynamic limit as long as the system remains fluid, i.e., provided no 
condensate appears, as could be expected. 

With the rate (jl.5[> . we have 

r(6 + i)fc! r 1 kl , , h , r(6 + i) , . 
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hence 



r 1 hl\ _ u )b-i 

P{z)= I -L ) —u = 2 F 1 (l,l;b+l;z), (3.5) 

1 — zu 



where 2 Pi is the hypergeometric function. 

The following values will be needed in the sequel: 



6-1' w (b- 1)(6- 2)' 

P "« - ( t -l)( t -2)( t -3) - <"> 

where it is understood that these quantities are convergent respectively for b > 1, 
6 > 2, and & > 3. More generally, the function P(z) has a branch cut at z = 1, with 
a singular part of the form 

(3.7) 

Whenever 6 is a positive integer or a half-integer, P(z) can be expressed in terms of 
elementary functions. In particular, for every integer b = n > 1, its singular part is of 

D sg( 

ln(l - z) 



the form P sg (z) ~ (— l)™n(l — z) n 1 ln(l — z). We have, for instance, 



6=1: P(z) 
6 = 2: P(z) 
6 = 3: P(z) 



z 

2(1 - z) ln(l - z) 2 

^2 



3(1 - z) 2 ln(l - z) 3 9 



*>(*)= K J r (3-8) 



z 3 z 2 2z' 

4(1 - z) 3 ln(l - z) 4 10 22 

z 3 z 2 3z 

The homogeneous ZRP with rate (|1.5[) has a continuous phase transition when 
the fugacity z reaches the singular point z c = 1. For 6 > 2, this takes place at a finite 
critical density (see (|1.6p ) 

* - 5$ " i^" (M) 

This critical density separates a fluid phase (p < p c ) and a condensed phase (p > p c ), 
whose main characteristics are as follows. 

• Fluid phase (p < p c ). The canonical and grand-canonical ensembles are equivalent 
throughout the fluid phase. The density p = p GC increases from to p c as the 
fugacity z increases from to 1. The occupation probabilities ff. = ff? c (see (|3.2[0 
fall off exponentially. 

• Critical density (p — p c ). This density is reached when the fugacity z takes the 
singular value z c = 1. It is therefore the maximal density that can be reached in 
the grand-canonical ensemble. 

At the critical density, the occupation probabilities 

_ Pk ^ (6-i)r(6) 
fk - p(T) ¥ — (3 - 10) 
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fall off as a power law with the continuously varying exponent b. The second 
moment of the occupation probabilities at criticality, 



is convergent for b > 3 (regime of normal fluctuations), whereas it is divergent 
for 2 < b < 3 (regime of anomalous fluctuations). 
• Condensed phase (p > p c ). The condensed phase only exists in the canonical 
ensemble where the total number N of particles is imposed. A large ZRP in its 
condensed phase consists of a uniform critical background, characterized by the 
critical occupation probabilities (13.101) . and of a single macroscopic condensate, 
containing an extensive number of excess particles of order A = M(p — p c ) 
(see (|1.3|l ). For a system where all sites are equivalent, the condensate can be at 
any site m with probability l/M. The fluctuations of the number of particles in 
the condensate around A are known to be Gaussian and to scale as M 1 / 2 in the 
regime of normal fluctuations (6 > 3), whereas they have a broad distribution and 
scale as M 1 ^ 6-1 ) in the regime of anomalous fluctuations (2 < b < 3) [TBI I15j . 

For the homogeneous ZRP, the spontaneous symmetry breaking (SSB) scenario 
in the condensed phase, with its unique condensate on top of a critical background in 
the thermodynamic limit, has been proved rigorously [131117] . The rare configurations 
where the condensate is not unique, i.e., where the A excess particles are shared by two 
sites in significant proportions, have been argued to occur with a probability falling 
off as the power law l/A/ 6-2 [14]. 

Let us close with a comment on universality. The condition b > 2 for the existence 
of a condensed phase is universal, in the sense that any rate with asymptotic behavior 
Uk = 1 + b/k + • • • yields a finite critical density p c only for b > 2. In contrast, the 
value of p c is not universal as it depends on the full shape of the rate Uk- Let us 
mention e.g. that the rate [51 [TUJ 



yields pk = (k + 1)~ 6 and p c = ((b — I) /((b), where £ is Riemann's zeta function. 
This choice of rate makes sense for all values of b, whereas (|1.5[) is limited to the range 
b > — 1 (as m = b + 1 > 0). The stationary state of the ZRP with this rate coincides 
with the static model analyzed in [7]. 

3.2. Occupation-independent inhomogeneous ZRP: explicit symmetry breaking 

We now consider the occupation- independent inhomogeneous ZRP (u^ — 1), where the 
single-particle weights q m are i.i.d. random variables drawn from the distribution (ll.4[) . 

The analysis of the model [3TJ 122] is again simpler in the grand-canonical 
ensemble. Taking the thermodynamic limit of (|2.32p . we readily obtain 



The critical density p c is then obtained by evaluating the above expression at the 
critical fugacity z c = l/q max = 1: 




(3.11) 




(3.12) 




(3.13) 




(3.14) 
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This integral is convergent for c > 1. It should however be underlined that the 
existence of a finite critical density p c is a characteristic of the thermodynamic limit. 
The density (|2.32p of any finite system indeed diverges, albeit at a further point, 
i.e., z — > 1/qi, where q± < 1 denotes the largest single-particle weight. For the 
distribution (|1.4[) we thus get (see (|l-7|l ) 

Pc = -^r. (3.15) 
c — 1 

The mean occupation distribution at the critical density is obtained by averaging 
the geometric law (|2.30[) over the distribution /(</), again at z = z c = 1. It reads 

t f 1 / x k ,, s cTic+l)k\ cT(c+l) 
We thus recover the property that the critical density 

Pc = J2 k fk ( 3 - 17 ) 

fc>0 

is convergent for c > 1. Finally, the mean square critical occupation 



r «(! + «) ^ c + 2 

>o Ji 

is convergent for c > 2. 

In the canonical ensemble, there is again a condensed phase for p > p c . In the 
latter phase, a large but finite system typically has a single condensate of size A 
(see (|1.3jl ). This condensate sits on the site with the largest single-particle weight (qi, 
say) with very high probability. This almost sure localization of the condensate on the 
most favored site is the gist of the explicit symmetry breaking (ESB) scenario at work 
in the occupation- independent inhomogeneous ZRP (see section B~2l for more details). 

The exponent c characterizing diffusion disorder therefore plays, in the present 
situation of ESB, a role analogous to that of the exponent b characterizing interactions 
in the situation of SSB on a homogeneous system. We have a fluid phase for c < 1 or 
b < 2, a condensed phase with anomalous fluctuations for 1 < c < 2 or 2 < b < 3, and 
a condensed phase with normal fluctuations for c > 2 or b > 3. 



4. Universal phase diagram 

We now address the study of the stationary state of the full inhomogeneous ZRP, 
characterized by the interaction exponent b entering the rate (|1.5[) and by the disorder 
exponent c entering the distribution (|1.4p of the single-particle weights. In this section, 
we aim at constructing the universal phase diagram of this ZRP in the plane of the 
exponents b and c, displayed in figures [1] and [2j 



4-1. Existence and nature of the condensate (critical density, critical fluctuations) 

The analysis is again easier within the grand-canonical approach. The mean critical 
occupation profile is obtained by taking the thermodynamic limit of (I2.16[) at the 
critical fugacity z c = 1. We thus obtain 
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and hence 



1 1 P '^ tt ^ 



^=l qP ' { % P " {q) ft«H (-) 

If the critical density p c is finite, the model has a condensed phase. On a 
large but finite system of M sites, if the particle density p = N/M exceeds p c , an 
extensive number A = M(p — p c ) of excess particles (see (|l-3|l ) will typically form a 
unique condensate. Furthermore the finiteness of p c will determine whether the size 
fluctuations of the condensate around A are normal or not. 

The uniqueness of the condensate, which has been proved rigorously in the 
homogeneous case [T3J E]i is basically due to the fact that the convexity of 
the interactions between particles suppresses the mixed configurations where the 
condensate is not unique, i.e., where the A excess particles are shared by two or more 
sites in significant proportions. This phenomenon is so general that it can already be 
exemplified on a ZRP on two sites for any b > (see figure IA1I in |Appcndix A[ ) . Let 
us mention that a rigorous combinatorial analysis of the condensation phenomenon in 
a finite system has been reported in [18] , 

It is therefore of central importance to control the convergence properties of the 
integral expressions (|4.2p for p c and p c . These properties are entirely dictated by the 
regime q — > 1. Two cases have to be dealt with separately. 

• For b > 1, -P(l) = b/(b — 1) is finite (see (|3.6|) ). As a consequence, the integral 
entering (|4.1j) falls off as k~ c . This observation yields the power-law decay 

f^^Tc ( b>1 )- (4-3) 

Thus p c is finite for b + c > 2, while p c is finite for b + c > 3. These two properties 
can be directly checked by analyzing the integrals (|4.2j) . They generalize the 
well-known conditions b > 2 and b > 3 for the homogeneous case, which are 
consistently recovered in the c — > limit. 

• For b < 1, P(q) ~ (1 — <z) b_1 is divergent (see (|3.7|) ). As a consequence, the 
integral entering (14.11) falls off as k ~ l . We thus obtain the power-law decay 

h-^T! (*<!)• ( 4 - 4 ) 

As a consequence, p c is finite for c > 1, while p c is finite for c > 2. These 
conditions, which were already derived in section[3]in the occupation-independent 
case (tife = 1), thus extend to the whole range b < 1. They can again be directly 
checked by analyzing the integrals (|4.2p . 

To sum up, the inhomogeneous ZRP defined by (jl.4[> and (jl.5l) exhibits a 
condensation transition only if attractive interactions and/or diffusion disorder excess 
a threshold value. The critical density p c is indeed finite for all c > if b > 2, for 
c>2 — &ifl<6<2, and for c > 1 if b < 1. Similarly, the critical mean square 
occupation p c is finite for all c > if b > 3, for c > 3 — b if 1 < b < 3, and for c > 2 if 
b < 1. The above inequalities demarcate the different regions of the phase diagrams 
of figures [T] and [2j 
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Whenever the critical density p c is finite, it is a rapidly decreasing function of the 
exponents 6 and c. The estimate 

p c « 1 (4.5) 

holds whenever 6 and c are simultaneously large. In the regime where c is small, for 
any b > 2, the critical density departs linearly from its value p c (0) = 1/(6 — 2) in the 
homogeneous case (see (|1.6p . (|3.9p ~). as 

(4.6) 



p c (c) =p c (0)(l-lc + -..), 



with 



1 = 



1 - 



gP'(g) 



/0 v p c (0)P(g)7 (4J) 

Figure [3] shows a plot of the ratio p c /p c (0) against c for 6 = 3 and 6 = 4. For these 
integer values of the exponent 6, the function P(z) admits the explicit forms Q3.8p . 
which ease the numerical evaluation of the integrals (I4.2[) and (|4.7[) . For 6 = 3, we 
have p c (0) = 1 and 1= 2.34704. For 6 = 4, we have p c (0) = 1/2 and X = 1.70735. 




0.2 0.4 0.6 0.! 



1.2 1.4 1.6 1.1 



Figure 3. Critical density p c , normalized by p c (0), against the diffusion 
disorder exponent c, for 6 = 3 and 6 = 4. 



4-2. Localization properties of the condensate 

It is intuitively clear that the site with the largest single-particle weight (qi , say) has 
the highest chance of hosting the condensate. However, if the spacings between the 
largest weights q m are small enough, it is plausible to expect a competition amongst 
the corresponding sites to host the condensate. The aim of this section is to provide 
a quantitative treatment of this issue. 

We first need to characterize the law of spacings between the q m . This can be 
achieved by using the order statistics of the original unordered variables. Let us 
reorder the sites according to decreasing single-particle weights q m and relabel them 
so as to have q\ > q2 > ■ . ■ > qu- A first estimation of the m-th ordered weight q m 
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can be obtained by equating the ratio m/M and the cumulative probability for the 
random variable q to be larger than q m , i.e., F(q m ), with 

F{q)= I /(<?') q' = (i-<z) c - (4.8) 



Jq 

We thus obtain the estimate 

A more precise prediction can be obtained as follows. In the regime of most interest 
where M is large, while the order m is kept finite, the set of rescaled variables 

x m = MF(q m ) = M(l - q m ) c (4.10) 

is asymptotically distributed as Poissonian points with unit density on the positive 
real line. This property is clear for c = 1, where the q m are uniformly distributed, 
so that F(q) = 1 — q. In this case, the variables obtained by ordering M 

points drawn independently and uniformly in the large interval [0,M]. The sequence 
thus obtained converges to Poissonian points with unit density. The general case then 
follows by changing variables from q to 1 — F(q). A full account of the connections 
between extreme-value statistics or order statistics and point processes can be found 
in [31 Sec. 7.3]. 

For Poissonian points with unit density, we have 

x m =Ti+T 2 ^ hr m , (4-11) 

where the distances r m = x m — x m ~\ between successive points are i.i.d. exponential 
variables with the distribution f T (T) = c~ T (and xq — 0). We thus obtain the more 
precise estimate 

(4.12, 

where the distribution of Gamma distribution T(to, 1): 

x m-l „— x 

= ^— ryr (4.i3) 

The mean and the variance of x m thus read 

X m = VBXXm — to- (4-14) 
So, as the order m gets large, the distribution of x m becomes peaked around m, with 



Gaussian fluctuations of relative order 1/y/m. In this regime f|4. 1 2[) simplifies to (14.9[) . 

We can now turn to the localization properties of the condensate in the stationary 
state of a large but finite system. This aspect of the problem is more subtle than the 
issued discussed in section 14.11 Indeed, a condensate can only be dealt with in the 
canonical ensemble at fixed AT, which is less amenable to analytical studies. Moreover, 
we are interested in properties of single sites (to = 1,2, and so on) and not in averaged 
properties over the whole system. We therefore have to rely on the following more 
heuristic line of reasoning. 

Consider a large but finite system in its condensed phase (M ^> 1, p > p c ), for 
b > and for a given draw of the single-particle weights q m . We already know the 
following: in a typical configuration, the condensate is unique, and its size fluctuations 
around A = M(p — p c ) are negligible. Our goal is to estimate the probability II m that 
site to hosts the condensate in the stationary state. The expression (|2.4p of the weight 
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of an arbitrary configuration shows that this probability (unnormalized so far) can be 
estimated by 

n m ~ qt (4-15) 
hence, using (|4.12|) . 

n m ^exp^A(^) 1/C ). (4.16) 

The argument of the exponential is proportional to A/M 1 ^ ~ M 1_1 / c . This estimate 
singles out the borderline value c — 1, where the dependence on the system size M 
drops out. This phenomenon already takes place in the simple case of the ZRP on 
two sites. In this situation, analyzed in | Appendix A| quantities of interest depend on 
the rescaled inhomogeneity parameter 9 = Ne (see (|A.5|) ). In the present case of a 
large system, the role of e is played by the difference between the largest two weights, 
i.e., qx - q 2 ~ M~ x l° (see (02J)), so that 9 ~ A/M 1 ^ ~ M 1_1 / c scales exactly as the 
argument of the exponential in (|4. 16[) . 

The following regimes have to be dealt with separately, according to the value of 
the exponent c with respect to the borderline value c = 1. 

• Extended regime (c < 1). In this first regime, p c is finite for b > 2 — c. The 
distribution f(q) of the single-particle weights diverges as q — !• q m ax = 1. Hence 
the few largest weights q\, q%, (?3, . . . pile up near q = 1. The typical distance 
between them, scaling as M _1 ' c , is much smaller than the mean spacing 1/M 
in the bulk of the distribution. As a consequence, the probability II m for site m 
to host the condensate takes appreciable values for many values of to. In the 
regime of interest, i.e., 1 <C m -C M, we have x m ps to, and so the estimate (|4.16l) 
simplifies to 

/ . / TO \ Vc 

n m ~ exp {-A (-) 

~ exp f-(p - p^M^^m 1 ^) . (4.17) 



The typical number R of favored sites can be estimated as the value of to such 
that the argument of the exponential in (|4.17l) is of order unity. We thus obtain 

M M 1 - 

R ~^ = l W 418 

A c [p - p c f 

In a typical configuration of the stationary state of the system, there is a unique 
condensate sitting at a well-defined site. The probabilities II m that the various 
sites to host the condensate are however extended over a large hosting set of 
favored sites. The typical size R of this set grows according to the sub-extensive 
law (I4.18|) . In the c — > limit we recover an extensive growth R ~ M, in 
agreement with the fully extended nature of the condensate in the homogeneous 
ZRP, with its spontaneous symmetry breaking mechanism. In the opposite regime 
(c —> 1), the result (I4.18|) crosses over to a localized condensate, in agreement with 
an explicit symmetry breaking mechanism. To sum up, the extended condensed 
phase, with its sub-extensive hosting set, realizes a continuous interpolation 
between the two condensation scenarios which were known so far. 
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• Borderline case (c = 1). In the borderline situation corresponding to the critical 
value c = 1 (and so p c is finite for b > 1), any dependence on the system size 
drops out of the estimate (I4.16|) . which reads 

n m ~e-^>™, (4.19) 

where the Poissonian points with unit density. The probabilities II m 

therefore typically take appreciable values on a small and fluctuating number of 
sites. | Appendix B| is devoted to a detailed study of these fluctuating probabilities. 
Distributions of this sort seem to have been first described in the context of 
randomly breaking an interval [35 . They have since then been met in several 
circumstances. A useful tool to investigate them is provided by the quantity Y 
(see (|B.3[) V which is somehow similar to the participation ratio used in the theory 
of Anderson localization [3B]- In the present situation, this quantity Y will be 
shown to have a non-trivial distribution, which is universal in the sense that it 
only depends on the parameter r = p — p c . 

• Localized regime (c > 1). For c > 1 (and so p c is now finite for all 6), the 
distribution f(q) goes to zero as q — > <7 max = 1- The few largest weights q±, q%, 
93, . . . are therefore rather distant from one another. The successive distances 
between them indeed scale as M -1 / c . As a consequence, for a typical draw of 
the single-particle weights, the first probability III is overwhelmingly larger than 
the other ones. The condensate is therefore localized on the most favored site 
(m = 1) with very high probability. The events where the largest two weights 
IIi an d n 2 become comparable occur with a small probability, which can be 
estimated to be of order M* 1— 1 / c ). These rare events will be analysed more 
precisely in |Appendix C| 

lf.,3. Phase diagram 

Putting all the above results together, we arrive at the universal phase diagram 
depicted in figures [I] and [2] The phase diagram shown in figure [T] exhibits three phases 
in the plane of the exponents b (characterizing interactions) and c (characterizing 
diffusion disorder). In the fluid phase, there is no condensate at any density. In the 
localized condensed phase, for p > p c , the condensate lives on the most favored site 
(m = 1) with very high probability. In the extended condensed phase, for p > p c , 
the condensate lives on one of the sites of a large set of favored sites, whose size R 
grows sub-extensively, according to (j4.18[) . Figure [5] shows a more detailed phase 
diagram, where the condensed phases are further separated into phases where p c is 
finite, hence the size fluctuations of the condensate around A are normal (i.e., Gaussian 
and growing as M 1 / 2 ), and phases where p c is divergent, hence those fluctuations are 
anomalous, i.e., not Gaussian and growing faster than M 1 / 2 . More precisely, these 
fluctuations scale as M 1 / for b < 1 and 1 < c < 2, and as AfVC'+c-i) for b > 1 
and 2 — b < c < 3 — 6. These growth laws can be respectively deduced from the 
expressions (|4.4[) and (|4.3|) . 

It should be clear from its derivation that the above phase diagram is universal, as 
announced in the introduction, in the sense that it holds for any rate whose asymptotic 
behavior at large k is of the form ~ 1 + b/k, hence P sg (z) ~ (1 — z) b_1 , and for any 
distribution of single-particle weights with a finite maximum g max and a power-law 
scaling as f(q) ~ (q max - q)^ 1 ■ 
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5. Numerical illustrations 

This section is devoted to numerical illustrations of the predictions made in the 
previous section. 

5.1. Density and occupation probability profiles 

Let us start at a rather qualitative level, and look at the growth of the density profile 
across a finite system as particles are added one by one. In order to avoid irregularities 
due to fluctuations, we take the deterministic ordered single-particle weights 

/2m - 1 \ 1/c 

9m = 1 -[ J hr) (m=l,...,M), (5.1) 

corresponding to the most uniform non-random sampling of the distribution (|1.4j) . 

The left-hand panels of figure |4] show plots of the canonical densities p rn = (N m ) 
of the successive sites (m = 1, 2, ...), against the mean density p — N/M of the 
system, for M — 50 and TV up to N max = 100, i.e., p max = 2, for a fixed interaction 
exponent 6 = 4, and the deterministic single-particle weights (|5.ip with three values 
of c. The densities p m are calculated by means of (|2.7p - (|2.9p . 

The ordering q\ > q 2 > qz > ■ ■ ■ of the weights always reflects itself into the same 
ordering for the densities, i.e., pi (black) > p 2 (red) > ps(green) > ... The various 
phases of the system are expected to manifest themselves as different asymptotic 
growth laws in various parts of the density profile. For c = 0.5, the densities of the 
first few sites are observed to grow at similar rates, at least at low enough density. 
For c = 2, only the most favored site has a rapidly growing density, whereas all the 
other densities remain microscopic. Finally, the situation for c = 1 looks intermediate. 
These findings are in agreement with the predicted transition at c = 1 between an 
extended and a localized condensed phase. As a consequence of the interactions, all the 
densities p m except the first one (pi) present a maximum for intermediate densities 
and eventually decrease at larger mean densities. This phenomenon is of the same 
nature as the overshoot observed in the impurity problem [24] . 

The right-hand panels of figure |4] show logarithmic plots of the local occupation 
probabilities f m< f. of the successive sites (to = 1, 2, . . .), against k, for the same model 
parameters with TV = 100 particles on M = 50 sites. These occupation probabilities 
are calculated by means of (|2.8j) . The favored sites, corresponding to the smaller 
values of m, present a non-monotonic profile of occupation probabilities, with a local 
minimum at k m i n ~ A/2, and a secondary peak at A: max w A, representing the 
condensate. Only the data for those favored sites have been plotted for clarity. The 
numbers of curves (20, 13, and 8) therefore give the numbers of those favored sites. 
For a given favored site to, the area under the secondary peak, 

N 

n m = J2 ( 5 - 2 ) 

fc=femi„ 

gives an operational estimate of the hosting probability of site to. As c increases, the 
number of favored sites decreases, and the hosting probabilities ±l m decrease more 
and more rapidly (remember the vertical scale on the right-hand panels of figure |4] is 
logarithmic) . 

The role of the disorder exponent c is further emphasized in figure [5j showing an 
enlargement of the first two datasets presented on the right-hand side of figure 01 in 
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Figure 4. Left: plots of the local densities p m = (N m ) at sites m = 1, 2, . . . 
against the mean density p = N/M, for M = 50 sites and up to 7V max = 100 
particles, with b — 4. Each panel corresponds to a value of c. Each color 
consistently corresponds to a site label (black for m = 1, red for m = 2, 
and so on). Right: logarithmic plots of the occupation probabilities / m ,fc 
of the favored sites (see text) against particle number k, for M = 50 and 
iV = 100. Other parameters and colors are the same as for the left-hand 
panels. 



the region of the condensate peak and on a linear scale. The contrast between the two 
situations corroborates our prediction. For c = 0.5 (in the extended phase), a large 
number of favored sites exhibit a visible secondary peak, and thus host the condensate 
with an appreciable probability. For c = 1 (at the borderline between the extended 
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and localized phases), only a small number of sites (4 on the plot) exhibit a visible 
secondary peak, and the areas under the peaks decrease very rapidly. In both cases 
the maxima of the peaks roughly coincide with the predictions k = A = N — Mp c , 
i.e., k = 87 for c = 0.5 (p c = 0.2555) and k = 92 for c = 1 (p c = 0.1641) (vertical 
dashed lines). 



c=0.5 c=l 




Figure 5. Enlargement of first two datasets presented on the right-hand 
side of figure 21 in the region of the condensate peak and on a linear scale. 
Vertical dashed lines: thermodynamic predictions for the mean condensate 
size A. 



5.2. A quantitative study of the hosting probabilities 

We now turn to the most novel prediction of this work, namely the existence of an 
extended condensed phase, where the condensate can be hosted by a large set of 
favored sites, whose size R obeys the sub-extensive growth law (|4.18|) . 

Let us first examine the individual hosting probabilities IT m . Considering c = 0.5 
for definiteness, the estimate (|4.17l) simplifies to the half-Gaussian law 

II m - exp I — I . (5.3) 

Figure [B] shows a logarithmic plot of the hosting probabilities, calculated according to 
the definition (|5.2I) . for all the sites m having a condensate peak, against m 2 /M, for 
b — 4 and p — 2, in three situations: a random sample of M = 100 sites, a random 
sample of M = 200 sites, and a deterministic sample of M = 200 sites with single- 
particle weights (|5.1j) . The black line with slope p — p c = 1.7445 shows the theoretical 
estimate (|5.3j) . The dataset of the deterministic sample is found to accurately follow 
the theoretical slope. The data for the two random samples exhibit the expected 
trend, although they are too noisy to be conclusive. This observation justifies our 
choice to use deterministic single-particle weights as illustrative examples in figures [4] 
and[SJ Finally, it is worth noticing that the sum of all the hosting probabilities, 

Si = ^n m , (5.4) 

n i 

is very close to unity in the three situations shown in figure [6] We indeed respectively 
find Si-1 = 7.92 x 10~ 6 , 4.73 x 10~ 6 , and 4.42 x 10~ 6 . The asymptotic sum rule 
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Figure 6. Logarithmic plot of the hosting probabilities Il m against m 2 /M, 
for b = 4, c = 1/2 and p — 2 in three cases (see text). Black line with slope 
p — p c — 1.7445: theoretical estimate (|5.3ll . 



We proceed by an investigation of the size R of the hosting set for the condensate. 
Along the lines of the theory of Anderson localization [36] , we can expect that a reliable 
estimate of R will be provided by the inverse participation ratio, 

i.e., the inverse of the quantity Y defined in (|B.3|) . This operational definition holds 
for any finite sample. The hosting probability H m of a favored site m is measured as 
the area under the condensate peak (see (|5.2[) ). and set to zero for the unfavored sites 
where the occupation probabilities f m ,k exhibit no such peak. The normalization in 
the numerator of (|5.5p is needed for the smaller systems, where the II m do not exactly 
sum up to unity. We now have a tool to help us revisiting the various regimes of the 
phase diagram. 

• Extended regime (c < 1). In this regime, it is legitimate to make use of 
the expression (|4.17|) for the hosting probabilities H m in order to estimate the 
quantity R defined in (|5.5|) . Evaluating the sums as integrals, we end up with the 
result that R is asymptotically self-averaging and grows as 



R^AM 1 -', A=[ r(c+l). (5.6) 

VP- Pc) 

The scaling estimate (|4.18p has thus been turned into a quantitative prediction. 
Figure [7] shows plots of the mean size R of the hosting set, for p = 2 and 6 = 3 and 
6 = 4, against M 1 / 2 for c = 0.5 (left) and against M 1 / 4 for c = 0.75 (right). All 
the datasets exhibit a linear behavior, corroborating the power-law growth ()5.6|) . 
and show a reasonable quantitative agreement with the predicted amplitudes A. 

Borderline case (c = 1). This borderline situation is investigated in detail 
in | Appendix B| The quantity R, which identifies with 1/Y (see (|B.3|) ). has a 
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c=0.5 c=0.75 




Figure 7. Mean size R of the hosting set in the extended regime. Straight 
lines have the theoretical slopes A (see (|5.6p ). Left: Data for c = 0.5 and 
p = 2 against M 1/2 . For 6 = 3, p c = 0.4140 and A = 0.9952. For b = 4 
(data shifted vertically for readability), p c = 0.2555 and A — 0.9489. Right: 
Data for c = 0.75 and p = 2 against M 1/4 . For 6 = 3, p c = 0.3078 and 
A = 1.0418. For b = 4 (data shifted vertically for readability), p c = 0.2007 
and A = 0.9949. 



non-trivial limit distribution which only depends on r = p — p c . An expression 
for the mean value R is given in (|B.31I) . Figure [8] shows a plot of R against M for 
p = 2 and b = 3 and 6 = 4. The data seem to exhibit a fast convergence to well- 
defined limits. The apparent limiting values however show a slight disagreement 
(of the order of half a percent in relative value) with the theoretical predictions 
given in (lB~3l"]) (horizontal lines), i.e., R = 1.9411 for b = 3 (p c = 0.2421, 
r = 1.7579), and R = 1.8995 for b = 4 (p c = 0.1641, r = 1.8359). The observed 
discrepancy is most certainly attributable to a very slow convergence. The data 
(not presented here) with the deterministic sampling (I5.1[) indeed exhibit a similar 
slow convergence, albeit toward different limiting values. 

Localized regime (c > 1). In this regime, we argued in section [4] that the 
condensate is localized on the most favored site (m = 1) with very high 
probability, hence R should converge to unity for large systems. The first 
correction to this trivial limit comes from the rare events where the weights III 
and II2 are comparable. The probability for this to occur falls off as 1/M 1_1 / C . 
In |Appcndix C| this argument is turned into the quantitative prediction (see (IC.7p ) 

fl-l, ^ M-NM. (5.7) 
2(P-Pc) 1 ; 

For c = 2 this yields 

- 7T 3 / 2 

iZ-l« =. (5.8) 

2{p- Pc )VM 

Figure [S] shows a plot of R against M for b = 4 and c = 2 (and so p c = 0.0927), 
and p — 2. The data exhibit a smooth maximum near R w 1.39 for M = 7, 
whereas they slowly go to unity for large systems, in good agreement with the 
theoretical prediction (f5T8|) . i.e., + 1A597/VM. 
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Figure 8. Mean size R of the hosting set in the borderline case (c = 1) 
against M, for p — 2 and fo = 3 and b — 4. Horizontal lines: predicted 
limiting values 1 = 1.9411 for b = 3 (p c = 0.2421, r = 1.7579), R = 1.8995 
for b = 4 (p c = 0.1641, r = 1.8359). 
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Figure 9. Mean size of the hosting set in the localized regime (c = 2) 
against M for b = 4 (and so p c = 0.0927), and p = 2. Full red line: 
theoretical prediction ifBTH)) . i.e., w 1 + 1.4597/VM. 



6. Complex networks 

In this section we investigate to what extent the above results concerning localization 
properties of the condensate apply to the ZRP on complex networks [37]. These 
networks are usually scale- free, i.e., characterized by a broad distribution of the node 
degrees K m , falling off as the power law Jk ~ K~ 7 , with 7 > 2. 

In the simplest situation where individual particles perform ordinary random 
walk, the stationary weight q m at node m is proportional to the degree K m . The 
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ZRP with this underlying inhomogeneous diffusion has been investigated by several 
groups [35J [2H1 [37] . Strictly speaking, this situation does not fit within the present 
work, as node degrees are unbounded. In other words, we have g max = oo and 
z c = 0, and so the model does not have a fluid grand-canonical phase. It has indeed 
been shown [25] that this inhomogeneous ZRP exhibits complete condensation (at 
least when Uk = 1), in the (weak) sense that the critical density p c vanishes in the 
thermodynamic limit. 

Consider a large but finite network consisting of M nodes. We have thus 
approximately A ps N = Mp. The nodes which can host the condensate are those 
with highest degrees, whose statistics has been studied recently [38] [39]. The node 
whose degree K\ is the highest is called the leader in those references. The estimate 
Ki ~ M x /(t-i) is easily obtained from an argument of extreme-value statistics. Two 
situations have to be dealt with separately. 

• There is a single leader, i.e., a single node with highest degree K\. In that case, 
we have K 2 = K\ — j for some j — 1,2,... and so —111(112/11!) rj jA/Ki ~ 
Mw~ 2 " vY- 1 ), The hosting probability il 2 is therefore exponentially small. The 
condensate is therefore hosted by the unique leader with an overwhelmingly high 
probability. 

• There are more than one node with highest degree K±. Let C = 2, 3, . . . denote 
the number of these co-leaders. In that case, the condensate is hosted by each 
co-leader with equal probability, and so we have R — C. It has been shown in [3j|] 
that the probability of having C co-leaders decreases rapidly with C, and that 
the leading event, namely the presence of two co- leaders (C — 2), occurs with a 
probability of the order of 1 jK\ . We thus obtain the estimate 

R- 1 ~ Nr 1 /^-^. (6.1) 

There is a striking similarity between this result and the prediction (|5.7[) in the 
localized regime (c > 1). The ZRP with inhomogeneous diffusion on complex networks 
is therefore always in the localized regime. The quantitative correspondence between 
exponents reads 

c=T- 2 , ,e., 7 =_. (6.2) 
The predictions made here are entirely due to the possible existence of co-leaders. 



7. Discussion 



In this work we have investigated the combined effects of interactions and diffusion 
disorder on the condensation phenomenon in the inhomogeneous ZRP. Our main 
findings have been summarized in the phase diagrams shown in figures [1] and [2] These 
universal phase diagrams are drawn in the plane of two exponents, the interaction 
exponent b of the rate (|1.5[) and the disorder exponent c of the distribution (|1.4[) of 
the single-particle weights modeling inhomogeneous diffusion. 

The most prominent feature put forward in this work is the existence of an 
extended condensed phase. This novel phase, which corresponds to the domain where 
< c < 1 and b > 2 — c, therefore shows up as a result of the combined effects of 
strong enough interaction and weak enough disorder. In this phase, a typical high- 
density configuration has a unique condensate on top of a critical background, but the 
condensate may be located at any site of a large hosting set of favored sites, whose 
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size grows sub-extensively as R ~ M l ~ c . The extended condensed phase can therefore 
be viewed as implementing a continuous interpolation between the two symmetry 
breaking scenarios known so far, namely SSB in the homogeneous ZRP and ESB in 
the occupation-independent inhomogeneous case. 

It is worth underlining the qualitative difference between the extended condensed 
phase of the inhomogeneous ZRP emphasized in this work and several alternative 
condensation scenarios in interacting particle systems which have been put forward in 
the recent literature. The homogeneous ZRP with non-conventional hopping rates Uk, 
exhibiting either a non-monotonic dependence on the occupation A: or a weak explicit 
dependence on the system size, may have several 'condensates' in its stationary 
state [40j [41]. These new scenarios include the possibility of having an extensive 
number of 'condensates' whose typical size is finite but parametrically large, or a sub- 
extensive number of 'meso-condensates', each of them having typically a sub-extensive 
population. Some mass transport models, where the hopping rate (generalizing Uk) 
also depends on the occupation of neighboring sites, have been shown to develop an 
extended condensate, i.e., a high-density structure containing all the excess particles 
and extending over a number of sites of order M 1 / 2 [42]. Another variant of the same 
class of models in one dimension leads to a phenomenon of explosive condensation, 
giving rise to a moving condensate sweeping the system at an accelerated pace [431 . 
The 'target process', introduced in [33], is dual to the ZRP, in the sense that the 
hopping rate now depends on the occupation of the arrival site. In the case of 
asymmetric dynamics in two dimensions and above, the stationary state again contains 
a one-dimensional extended condensate, aligned with the direction of the mean current. 
The 'inclusion process', where particles are allowed to make non-local jumps over 
the lattice, exhibits unusual types of complete condensation, including a case where 
almost all the particles are condensed on the right-most side of a finite one-dimensional 
chain 45i . Finally, a non-Markovian generalization of the ZRP, where the motion of 
particles is governed by internal clocks endowed with their own dynamics, leads to 
a variety of non-equilibrium stationary states. In the asymmetric one-dimensional 
situation, the model may exhibit a localized condensate sitting on two neighboring 
sites and moving ballistically [4"5] . 

The present work focussed on the stationary state of the ZRP with inhomogeneous 
diffusion (|1.4[) and interaction (|1.5j) . The dynamical consequences of our findings are 
clearly the next question of interest. For the case of homogeneous diffusion on a large 
but finite system, the ergodic motion of the condensate in the stationary state has 
been investigated by the present authors in 14 . The leading mechanism turns out 
to be that all the excess particles sequentially quit the condensate and progressively 
rebuild it at another random distant site. The associated characteristic time scale can 
therefore be analyzed by means of an effective two-site model, with the key ingredient 
being the occupation probability profile fk of the latter model. This ergodic time scale 
is found to grow with the system size as the power law 

Ab+l 

r ~ ~ M\ (7.1) 

i.e., faster than the diffusive scale M 2 in the thermodynamic limit (as b > 2), but not 
exponentially fast. (In the special situation of symmetric dynamics in one dimension, 
our prediction is increased tor ~ A &+1 .) The scaling prediction (|7.ip has been alluded 
to or used in several situations germane to the present one [25] [301 EH 1331 1331 133 HS] ■ 
It has also been corroborated by a recent rigorous analysis [191 . albeit in the 
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regime where N becomes large while M is kept finite. Coming back to the 
inhomogeneous ZRP studied in the present work, the occupation probability profile on 
two sites exhibits a universal 'dip', scaling as N~ b , irrespective of the inhomogeneity 
(see (|A.17|1 ). It can therefore be expected that the scaling law (|7.1|) will still hold all 
over the extended condensed phase of the model (for b > 2), where the condensate 
is allowed to live on a large hosting set of favored sites. In the region of the phase 
diagram where 2 — c < b < 2 the ergodic time should scale as the diffusive one. Finally 
the coarsening process of the formation of the condensate is yet another facet of the 
dynamics worth investigating. 
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Appendix A. The case of two sites 

Let us consider the simplest non-trivial situation of the inhomogeneous ZRP on two 
sites. The quantities of interest can be given explicit expressions in this situation, while 
its behavior is already illustrative of many features of the model on larger systems. 

The single-particle weights qi and <?2 and the hopping rates W12 and W21 are 
related by W12Q1 — (see (|1.1[> ). We introduce an inhomogeneity parameter e so 

that q\jq2 — W21/W12 = e 2e , and parametrize the single-particle weights as 

qi = e e , <?2 = e~ e . (A.l) 

For definiteness, we set e > 0, so that q\ > q2, i.e., site 1 is favored. 



Appendix A.l. The occupation-independent case (u^ = 1) 

Let us start with the occupation-independent case (v,k = 1), and consider the canonical 
ensemble where the total number N of particles is fixed. 
The partition function (|2.36p becomes 

- g^ +1 Binh(iV + l)e 

^2,N = = T-j ■ \ A - Z ) 

qi — q2 sinhe 
Setting Ni — k and N2 = N — k, the occupation probabilities read 

aka N ~ k sinhe e ( 2k - N ) s 

fk = Prob^i = k} = -^r = , ' n ■ (A.3) 

Z2, n smh(A< + 1)6 

The occupation probability profile is thus exponentially increasing from the least 
probable event (site 1 is empty, i.e., k = 0) to the most probable one (site 2 is empty, 
i.e., k = N). The ratio between the smallest and the largest of these probabilities 
reads 



N 

_ „-2Ne 



(A.4) 



fa _ ( gg 
In \qi. 

When the number of particles N becomes large, in order for the populations of 
the two sites to remain comparable, one has to keep the model in the regime of a weak 
inhomogeneity, where e scales as 1/N. In this regime, many results assume scaling 
forms, involving the rescaled inhomogeneity parameter 

6 = Ne. (A.5) 
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For instance, the normalized occupation probability profile becomes 

(2x-l)0 

fk « sinhg (0 < x = k/N < 1) (A.6) 

and, in particular, the probability ratio (|A.4|) becomes e~ 2e . The density contrast 

d 

5p = Pi-P2 = (Ni -N 2 ) = — In Z 2 , N (A.7) 

scales as 

6pnN9> {6), (A.8) 
where the scaling function 

$ (6>) = cotanhfl - \ (A.9) 

starts increasing linearly as &o(0) ~ 0/3 for small 0, and slowly saturates to unity as 
$o(0) » 1- 1/0 at large 0. 

Let us briefly consider the grand-canonical case. We have (see (|2.3ip ) 

p GC = ^_ p GC = ^L_ (A . 10) 

f — ze e 1 — ze E 

The fugacity z is determined by imposing the mean total occupation (N GC ) = 
p GC +p GG . Skipping details, we again find that the density contrast 5p GC = p GC — p GC 
obeys a scaling law of the form (|A.9|) . with 9 = (N GC )e and 

$? C (0) = j . (ATI) 

This scaling function is rather close to its canonical counterpart (|A.9j) . It starts 
increasing as & GC (9) 0/2 and slowly saturates to unity as & GC (9) » 1 — 1/9. Both 
functions will be plotted in figure |X2j 

Appendix A. 2. The general case 

We now turn to the general case of the inhomogeneous ZRP on two sites with the 
rate (|1.5j) . considering only the canonical ensemble and keeping the same notations as 
above. 

The occupation probabilities read 

fk = Prob{Ai = k} = PkPN - k , (A.12) 

^2, N 

where the partition function is 

N 

z 2 , N = j2 ei2k ~ N)£ p k P N - k - ( A - 13 ) 

fc=0 

In the presence of attractive interactions, i.e., for b > 0, the most probable event 
is still that site 2 is empty (i.e., k — N), and the result (|A.4j) still holds. More 
interestingly, however, the occupation probability profile approximately reads 

„, r(6 + iy c("-")« 

/fc ~ Z 2 ,n k\N-k) b ' [ ' 
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whenever both populations are large (i.e., fc> 1 and N—k 3> 1). The above expression 
is the product of an increasing exponential (reflecting the inhomogeneity) and of two 
decaying power laws (reflecting the attractive interactions). 

In the scaling regime of a weak inhomogeneity (N large, e small, 8 = Ne fixed), 
the above probability profile has a deep minimum for a non-trivial value fc m j n of the 
occupation k, namely 

N 

kmin~A(e)-, (A.15) 



with 



A{8) 



2b 



b + 8 + y/b 2 + 8 2 
The corresponding smallest probability scales as 

/mill 



N 



B(8) 2 T(b + 1) N~ 



with 



B(0) 



b + y/b 2 + 8 2 
2b 



(A.16) 



(A.17) 



(A.18) 



The probability profile therefore exhibits a universal dip, in the language of [TJ] , whose 
scaling in N~ b holds irrespectively of the rescaled inhomogeneity parameter 8, which 
only enters the scaling functions A(8) and B(9). These functions are normalized so 
that A(0) — -B(O) = 1 in the homogeneous limit. The occurrence of a universal dip in 
the occupation probability profile is the gist of the uniqueness of the condensate in the 
inhomogeneous ZRP for any b > 0. It also has far-reaching dynamical consequences, 
which are briefly discussed in section [7] Figure IA1I shows a logarithmic plot of the 
occupation probabilities fk for N = 100 particles on two sites, for 6 = 1 and several 
values of 8. 
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Figure Al. Logarithmic plot of the occupation probabilities fk for iV = 100 
particles on two sites, with 6=1 and several values of 9. 



In order to derive more quantitative predictions concerning the scaling behavior 
of the partition function ^2,iv an d of the density contrast Sp, the following two regimes 
have to be considered separately, along the lines of [14] • 
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For 6 > 1, the sum in (|A.13|) is dominated by configurations such that almost all 
the particles sit on one site, i.e., either k is finite, or N — k is finite. Indeed the 
sum -P(l) of the Pk is convergent (see (|3.6|l ). In other words, there is complete 
condensation in the N — > oo limit. 
We thus obtain in the scaling regime 

Z ^ * 2bT b ib - + l 1} N ~" C ° Sh9 (A - 19) 

and 

5pwJV$i(fl), (A.20) 

where 

$i(0)=tanh0 (A.21) 
is independent of b in the range b > 1. The above scaling function starts linearly 
as $i(0) ~ and saturates exponentially fast to unity as $i(0) «1- 2e . 
For b < 1, the sum of the pk is divergent, and so all values of k contribute to 
the sum in (IA.13[) . Evaluating this sum in the scaling regime as an integral over 
y = 2x — 1 = 2k/N — 1, we obtain 

Z 2 N » AT1-2& F (Q) (A.22) 

* r(2-26) sin 2 (^6) K h K ' 

with 

T(3/2-b) e^y 

bl j Vir(i (i -iff 

= T(3/2 - 6)(0/2) b - 1 / 2 / 1/2 _ 6 (0), (A.23) 
where I v {9) denotes the modified Bessel function of index v. The above scaling 
function is again normalized so that -Fb(O) = 1 in the homogeneous limit. We 
thus have 

5p&N$ b (d), (A.24) 

with 

This scaling function starts linearly as $f,(0) 0/(3 — 2b) and slowly saturates to 
unity as $ 6 (0) w 1 - (1 - 6)/0. The result (|A.9|) is recovered for 6 = 0. 
The function $6(0) coincides with the function describing Langevin paramag- 
netism in dimension d = 3 — 2b. Indeed, setting y = cos u, the integral (1A.23[) can 
be recast as 

F b (0) = ■= ^. ( y 2) . . . r(sinu) d -V cos "u, (A.26) 



/?rr((d- l)/2) 

where we recognize the normalized partition function for a classical spin, 
represented by a vector s describing the unit sphere in d dimensions, coupled 
to a reduced magnetic field 0, so that u is the angle between s and the direction 
of the field. This analogy makes sense as long as s has continuous degrees of 
freedom, namely for d > 1, i.e., precisely b < 1. 

Figure IA2I shows the scaling function of the density contrast for several b. 
This function is independent of b and saturates exponentially fast for b > 1, whereas 
it depends continuously on b and saturates slowly for b < 1. Furthermore, at least 
for the case 6 = 0, the grand-canonical scaling function is very close to its canonical 
counterpart. 
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Figure A2. Finite-size scaling functions $(,(#) describing the density 
contrast of the two-site problem, against the rescaled inhomogeneity 
parameter 8, for various values of the exponent b. The grand-canonical 
scaling function <E>p c (#) is also plotted for comparison. 



Appendix B. Universal fluctuations in the borderline situation (c = 1) 

This Appendix is devoted to a detailed study of the fluctuations of the hosting 
probabilities II m of a large system in the borderline situation (c = 1). We are 
thus led to consider an infinite ordered sequence of unnormalized weights of the form 
(see dUg) ) 



n, 



(m = l,2,...), 



(B.l) 



where the x m are Poissonian points with unit density, whereas r = p — p c is the control 
parameter of the problem. 

We are interested in the joint distribution of the two sums 



and chiefly in the distribution of the ratio 
S 2 



Y 



Sf 



(B.2) 



(B.3) 



introduced in [35J. 

The formula (|4. 1 1|) allows us to recast the definition of the sums Si and S2 into 
the form 

51 = Xi + X1X2 + X1X2X3 + ■ • • , 

5 2 = X\ + X\X\ + X\X 2 2 X\ + • • • , (B.4) 
where the X m = cr rTm are i.i.d. random variables with the distribution 

f x (X) = -X- 1+1 l r (0<X<1). (B.5) 
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Infinite random sums such as (|B.4I) are known as Kesten variables [4"Tl |4"5] . They have 
been shown to play a role in various areas of the physics of one-dimensional disordered 
systems [451 [SHI HI] • I n the present situation, we have the identities 

Si = X{l + S{), S 2 = X 2 {l + S' 2 ), (B.6) 

where X is drawn from the distribution (|B.5j) and S[ and S' 2 are copies of the variables 
Si and S 2 , independent of X. The above identities can be viewed as a stochastic 
dynamical system for two degrees of freedom, described by the coupled variables S± 
and S'2, and submitted to the same noise X. 



Appendix B.l. Distribution of the normalization sum Si 

The distribution of a Kesten variable such as Si can only be worked out for special 
distributions of the X m |50] . The present situation, where the X rn have a power- law 
distribution on the interval [0, 1], is the simplest of all the exactly solvable cases [48]. 
It is convenient to use the Laplace transform of the distribution, 



L 1 (i) = e- tS K (B.7) 
The identity (|R6)l for Si yields 

Li(t) = e-'^+s;) = f 1 f x (X)e- tx Li(tX)X. (B.8) 
Jo 

The change of variable from X to u — tX leads to 

1 r* 

Li{t) = -t- 1/r u - 1+1/r c~ u Li(u)u, (B.9) 
r Jo 

hence the differential equation 



L[(t) = — Li(t) (B.10) 

rt 

and the explicit expression 

L 1 (t) = e- Fl(t)/r , (B.ll) 

with 



Fi(t)= [\l 
Jo 



°) 

X 



1/0 n>l 

where 7 = 0.577 215 . . . denotes Euler's constant, and Ei is the exponential integral. 

The series in the rightmost side of (|B.12[) demonstrates that the cumulants of Si 
have the simple expression 

5f (c) = — . (B.13) 
nr 

Wc have in particular Si = \/r and varSi = l/(2r), and so 1/r can be interpreted as 
the effective number of significant weights which contribute to the sums (|B.2j) . 
The full distribution of Si reads 

fsASi)= [ -Le t *-*W/'\ (B.14) 

./ 27T1 
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Figure Bl. Distribution fs 1 (Si) of the normalization sum Si, for three 
characteristic values of the control parameter r. 



Figure IBll shows a plot of this distribution for three characteristic values of the 
control parameter r. Each dataset consisting of 200 bins has been obtained by means 
of 10 9 iterations of the random recursion (IB. 61) . 

The distribution fg 1 becomes a narrow Gaussian in the r — ► limit, in qualitative 
agreement with the law of large numbers, as the number 1/r of significant weights is 
large. It gets progressively broader as r increases. 

The behavior of fs 1 for small values of Si can be estimated by using large positive 
values of t in the contour integral (|B.14[) . Neglecting the exponentially small function 
Ei(— t), we readily obtain the power law 

fsASl)x mj^ sll+1/r &^°)- (B.i5) 

The density therefore diverges as Si — ► for r > 1, whereas it tends to zero for r < 1. 
This feature is clearly visible on figure |B~TI fg t vanishes linearly for r = 1/2, goes to 
the constant e~ 7 for r = 1, and diverge as l/y/Si for r = 2. 

The behaviour of fs 1 {Si) for large values of Si can be estimated from the contour 
integral (|B.14|) by means of the saddle-point method. Skipping details, we are left with 
a super-exponential tail of the form 

f Sl (Si) - cxp(-Si(lnSi + ln(lnSi) + lnr + •••)) • (B.16) 

Finally, performing an integration by parts in (|B.14[) . we obtain the following 
differential-difference equation for the distribution fs 1 (Si): 

(r- l)f Sl (Si) + rSif Sl (Si) + f Sl (Si - 1) = 0. (B.17) 

This functional equation has many consequences, including that the estimate (|B.15|) 
holds identically over the interval < Si < 1, and that the distribution /s^Si) has 
weaker and weaker singularities at all the integer values Si = 1,2,... Distributions 
of this kind have been met in several instances, starting with Renyi's analysis of the 
one-dimensional car parking problem |52j . 
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Appendix B.2. Joint distribution of the sums S\ and S 2 

The joint distribution fs 1 ,s 2 (Sii S 2 ) of the sums Si and S 2 can be studied along the 
same lines. Its Laplace transform reads 



L 2 (h,t 2 ) =e-*iSi-*2S 2 . ( R18 ) 
The identities (|B~6| yield 



L 2 (ti,t 2 )= [ fx(X)e~ tlX ~ t2x2 L 2 (tiX,t 2 X 2 )X. (B.19) 
Jo 

Setting g = t 2 /t 2 , the change of variable from Itoii = t\X leads to 

L 2 {t u gt\) = - ti 1/r f 1 u~ 1+1/r e~ u ~ 9u2 L 2 (u, gu 2 ) u. (B.20) 
r Jo 

This is an integral equation for the function L g (ti) = L 2 (ti, gt 2 ), which therefore 
obeys the differential equation 

L' g (h) = — L g (h). (B.21) 

We thus obtain the explicit expression 

L 2 {t 1 M)=e- F2(tlM),r i (B.22) 



with 

F 2 (h,t 2 ) = I (1-e"' 1 



1 

t\x— t^x 1 



( 1 \m-fn-Um+ii 

= y h h , (b.23) 

(m,n)^(0,0) v ; 

The above result has the striking feature that the generating function F 2 (ti,t 2 ) obeys 
the inverted heat equation 

dF 9 d 2 F 2 

!»7 + T«f = - < R24 > 

where t 2 plays the role of a negative time, and with the initial condition ^2(^1,0) = 
F\(ti). The series in the rightmost side of (|B.23j) demonstrates that the joint 
cumulants of Si and S 2 have the simple expression 

L 



STSf = 7 =-r-. (B.25) 

1 2 (m + 2n)r 



The joint distribution of Si and S 2 , 



f Sl MSl,S 2 )= I / iL^^Si^-F^M/r (B26) 
' ' 2ni 2ni 



is however highly non-trivial. 
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Appendix B. 3. Distribution of the quantity Y 

The distribution of the quantity Y introduced in (|B.3[) reads formally 

poo 

f Y (Y)= fsuS^SuY^SfS^ (B.27) 
Jo 

where the joint distribution fsi,S 2 (Sx, S2) of the sums Si and S2 is given by (|B.26[) . 
The resulting triple integral formula is not very useful. Figure IB2I shows a plot of the 
distribution fy(Y) for four characteristic values of the control parameter r. This 
distribution inherits from the joint distribution fsx.Sz (^l; $2) & highly non-trivial 
structure, with weaker and weaker singularities at the reciprocal integers Y = 1/2, 

Y = 1/3, and so on. As the control parameter r gets larger and larger, the overall 
distribution shifts towards larger values of Y, and its singularities at Y = 1 and 

Y = 1/2 become more and more pronounced. 




0.2 0.4 0.6 0.8 1 

Y 



Figure B2. Distribution fy{Y) of the quantity Y, for four characteristic 
values of the control parameter r. 



The moments Y k can be evaluated, at least in principle, by means of the identity 

ok ok poo 

Y k = %= - ^ - / e~* Sl t 2k ~ x t. (B.28) 
Sf k (2k -1)! J ■ K 1 

The quantities S k c~ tSl can be calculated by expanding the result (|B.22I) as a Taylor 
series in ti. Some algebra involving integrations by parts leads to the following integral 
expressions for the first two moments of Y: 

1 r°° 

Y = 1 - - c-^/'-e-'t, 
r Jo 

7* = 1 - — / e ~ F ^' r e-*(i 2 + 3t + 5)t 
6r In 



e 



- Fl W/ r e- i (2-(t + 2)e-*)t. (B.29) 



1 

6r 2 

Figure |B3] shows plots of Y (left) and of varF (right), obtained by evaluating the 
above integrals numerically, against r/(r + 1). The variance reaches its maximum 
vary « 0.05218 for r « 2.25. 
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Figure B3. Plots of Y (left) and of var Y (right), against r/(r + 1). 



The mean value of the quantity R = 1/Y can be determined along the same line 
of thought. Some algebra using the identity 

R=l = ^=S!j o e-^t (B.30) 

leads to the following integral expression 

poo i 

R=1 + ^L e - Fl W/W(erf^) 2 i. (B.31) 
4H J t 

In the r < 1 regime, where the number of significant weights gets large, we obtain 

— r r 2 r 3 101r 4 



2 12 24 720 



r 

var Y = — 



2 360 90 
2 1 14r 

r 3 15 

and so the distribution of Y is asymptotically a narrow Gaussian around r/2, whose 
width scales as r 3 / 2 . 

In the opposite r > 1 regime, the distribution of Y is dominated by a strong 
divergence as Y — > 1, and so the moments approach unity as 

— 1 In 2 

Y =1 -- + — + •••> 

^ 1 41n2 + 7 

varF = — +■■■, 

3?' 12r z 

— 7T 7rln2-2C 

R =l + _ + _ S3 _ + ... 1 (B.33) 

where C = 0.915 965 . . . denotes Catalan's constant. This regime is studied in more 
detail in the next section. 

Appendix B.4- Effective two-site model in the r» 1 regime 

In the r> 1 regime, the (ordered) weights Ii m — e~ TXm fall off very rapidly with m. 
We are therefore led to consider an effective two-site model, where only the largest two 
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weights are kept, while all the other ones are neglected. The corresponding Poissonian 
points are X\ = T\ and x 2 = T\ + t 2 . 
Within this framework, we have 

Si = n 1 + n 2 = e- rTi (i + c- rT2 ), 

S* 2 = II 2 + n 2 = e _2m (1 + e- 2rT2 ), (B.34) 

and hence 

1 + e~ 2r7 " 2 

Y = T —^ -. (B.35) 

The quantity Y therefore only depends on r 2 . It lives in the interval 1/2 < Y < 1, 
and its distribution reads 

1 / 1-Y V /r 

f Y {Y) = - — == . (B.36) 

1 v ' r(l - y)V2F - 1 \y + V2F - 1 / v ; 

This distribution exhibits an inverse-square-root singularity at the lower edge: 

M Y )*^Tf ^^V2), (B.37) 
and a very strong divergence at the upper edge: 

f Y (Y)^-(l-Y)- 1+1/r (Y->1). (B.38) 
r 

As a consequence, all the moments of Y approach unity according to 

yT=l_^ ; (B.39) 
r 

where the amplitude of the negative 1 jr correction term reads 
A y 1 - Y k 

The generating series of these numbers can be evaluated as 

Y 



k>i Jl 



1 - z J 1/2 1(1 - zY) 



Z (B.41) 



1 - z V 2 - z 1 - z 

This expression yields the values /zi = 1 and /i 2 = 5/3 (in agreement with (|B.33[) to 
order 1/r), ^3 = 32/15, (14 = 52/21, ^5 = 863/315, and so on, as well as the linear 
recursion formula 

(2k - 1H - (3fe - 2) Mfc _j + {k- l)/i fe _ 2 = 1, (B.42) 

and the asymptotic growth law fi/. ~ ln(2fc) + 7. 

Finally, the formal expression obtained by setting k = — 1 in (|B.40|) . i.e., 

agrees with the expansion (|B.33|) for i?. 
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Appendix C. Effective two-site model in the localized regime (c > 1) 

The effective model studied in |Appendix B| can be extended to the localized regime 
(c > 1), in order to describe the rare events where the first two weights 111 and II2 are 
comparable. We still have S\ = H\ + II 2 and S2 = flf + flf , and so 

1 + v 2 (1 + v) 2 2~n 

Y = - ^, R=\ ^ = 1 + ^ — L a, C.l 

(I+77) 2 I + 77 2 I + 77 2 v ' 

where the weight ratio 77 = ILi/IIi reads (see (14.161) ) 

^^(-^(^ -**/•)). (C.2) 
Let us first evaluate the mean value of this ratio: 

* * r ^ r exp (-^ - xi/c )) *• (c - 3) 

The integral over x\ is dominated by small values of the difference T2 = X2 — x\ ~ 

Linearizing the argument of the exponential with respect 
to T2, and integrating over T2 first, we readily obtain 

* « cr(2 - 1/c) ^ = ^^-^ M-(-V C ). (C.4) 

Higher moments of 77 can be estimated by replacing in the above result A by the 
product fcA. We thus get immediately 

(C.5) 

This result leads to the picture that 77 is zero with very high probability, and of order 
unity with a small probability of order 77. The quantity fj thus provides an operational 
estimate for the probability that ITi and II2 are comparable. This probability was 
already anticipated in section|4]to fall off as M~( 1_1 / c ). 

Finally, by expanding the expression (|C.1|) for R as a power series in 77, using (|C.5|) . 
as well as the resullQ 

we are left with the following quantitative estimate for the mean value of R in the 
localized regime: 

r - 1 « ™r(2-i/ c) *± = ^r(2-i/ c) (1 _ 1/e) 

2 A 2(p- Pc ) y ' 
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